*
* $Id$
*
* $Log: samcst.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:37  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:16  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:04  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
*-- Author :
*$ CREATE SAMCST.FOR
*COPY SAMCST
*
*=== samcst ===========================================================*
*
C
C*******************************************************************
C
      SUBROUTINE SAMCST(KPROJ,EKIN,CST)
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
C
C  HJM 24/10/88
C
C                 SAMPLING OF COS(THETA)
C                 FOR NUCLEON-PROTON ELASTIC SCATTERING
C                 ACCORDING TO HETKFA2/BERTINI PARAMETRIZATION
C
C-------------------------------------------------------------------
C
      DIMENSION DCLIN(195),DCHN(143),DCHNA(36),DCHNB(60)
      DIMENSION PDCI(60),PDCH(55)
      REAL RNDM(4)
      SAVE DCLIN, DCHN, DCHNA, DCHNB, PDCI, PDCH
C
      DATA (DCLIN(I),I=1,80) /
C***     DCLN ARRAY
     *     5.000D-01,  1.000D+00,  0.000D+00,  1.000D+00,  0.000D+00,
     *     4.993D-01,  9.881D-01,  5.963D-02,  9.851D-01,  5.945D-02,
     *     4.936D-01,  8.955D-01,  5.224D-01,  8.727D-01,  5.091D-01,
     *     4.889D-01,  8.228D-01,  8.859D-01,  7.871D-01,  8.518D-01,
     *     4.874D-01,  7.580D-01,  1.210D+00,  7.207D-01,  1.117D+00,
     *     4.912D-01,  6.969D-01,  1.516D+00,  6.728D-01,  1.309D+00,
     *     5.075D-01,  6.471D-01,  1.765D+00,  6.667D-01,  1.333D+00,
     *     5.383D-01,  6.054D-01,  1.973D+00,  7.059D-01,  1.176D+00,
     *     5.397D-01,  5.990D-01,  2.005D+00,  7.023D-01,  1.191D+00,
     *     5.336D-01,  6.083D-01,  1.958D+00,  6.959D-01,  1.216D+00,
     *     5.317D-01,  6.075D-01,  1.962D+00,  6.897D-01,  1.241D+00,
     *     5.300D-01,  6.016D-01,  1.992D+00,  6.786D-01,  1.286D+00,
     *     5.281D-01,  6.063D-01,  1.969D+00,  6.786D-01,  1.286D+00,
     *     5.280D-01,  5.960D-01,  2.020D+00,  6.667D-01,  1.333D+00,
     *     5.273D-01,  5.920D-01,  2.040D+00,  6.604D-01,  1.358D+00,
     *     5.273D-01,  5.862D-01,  2.069D+00,  6.538D-01,  1.385D+00/
      DATA (DCLIN(I),I=81,160) /
C***     DCIN ARRAY
     *     5.223D-01,  5.980D-01,  2.814D+00,  6.538D-01,  1.385D+00,
     *     5.202D-01,  5.969D-01,  2.822D+00,  6.471D-01,  1.412D+00,
     *     5.183D-01,  5.881D-01,  2.883D+00,  6.327D-01,  1.469D+00,
     *     5.159D-01,  5.866D-01,  2.894D+00,  6.250D-01,  1.500D+00,
     *     5.133D-01,  5.850D-01,  2.905D+00,  6.170D-01,  1.532D+00,
     *     5.106D-01,  5.833D-01,  2.917D+00,  6.087D-01,  1.565D+00,
     *     5.084D-01,  5.801D-01,  2.939D+00,  6.000D-01,  1.600D+00,
     *     5.063D-01,  5.763D-01,  2.966D+00,  5.909D-01,  1.636D+00,
     *     5.036D-01,  5.730D-01,  2.989D+00,  5.814D-01,  1.674D+00,
     *     5.014D-01,  5.683D-01,  3.022D+00,  5.714D-01,  1.714D+00,
     *     4.986D-01,  5.641D-01,  3.051D+00,  5.610D-01,  1.756D+00,
     *     4.964D-01,  5.580D-01,  3.094D+00,  5.500D-01,  1.800D+00,
     *     4.936D-01,  5.573D-01,  3.099D+00,  5.431D-01,  1.827D+00,
     *     4.909D-01,  5.509D-01,  3.144D+00,  5.313D-01,  1.875D+00,
     *     4.885D-01,  5.512D-01,  3.142D+00,  5.263D-01,  1.895D+00,
     *     4.857D-01,  5.437D-01,  3.194D+00,  5.135D-01,  1.946D+00/
      DATA (DCLIN(I),I=161,195) /
     *     4.830D-01,  5.353D-01,  3.253D+00,  5.000D-01,  2.000D+00,
     *     4.801D-01,  5.323D-01,  3.274D+00,  4.915D-01,  2.034D+00,
     *     4.770D-01,  5.228D-01,  3.341D+00,  4.767D-01,  2.093D+00,
     *     4.738D-01,  5.156D-01,  3.391D+00,  4.643D-01,  2.143D+00,
     *     4.701D-01,  5.010D-01,  3.493D+00,  4.444D-01,  2.222D+00,
     *     4.672D-01,  4.990D-01,  3.507D+00,  4.375D-01,  2.250D+00,
     *     4.634D-01,  4.856D-01,  3.601D+00,  4.194D-01,  2.323D+00/
C
      DATA PDCI /
     *     4.400D+02,  1.896D-01,  1.931D-01,  1.982D-01,  1.015D-01,
     *     1.029D-01,  4.180D-02,  4.228D-02,  4.282D-02,  4.350D-02,
     *     2.204D-02,  2.236D-02,  5.900D+02,  1.433D-01,  1.555D-01,
     *     1.774D-01,  1.000D-01,  1.128D-01,  5.132D-02,  5.600D-02,
     *     6.158D-02,  6.796D-02,  3.660D-02,  3.820D-02,  6.500D+02,
     *     1.192D-01,  1.334D-01,  1.620D-01,  9.527D-02,  1.141D-01,
     *     5.283D-02,  5.952D-02,  6.765D-02,  7.878D-02,  4.796D-02,
     *     6.957D-02,  8.000D+02,  4.872D-02,  6.694D-02,  1.152D-01,
     *     9.348D-02,  1.368D-01,  6.912D-02,  7.953D-02,  9.577D-02,
     *     1.222D-01,  7.755D-02,  9.525D-02,  1.000D+03,  3.997D-02,
     *     5.456D-02,  9.804D-02,  8.084D-02,  1.208D-01,  6.520D-02,
     *     8.233D-02,  1.084D-01,  1.474D-01,  9.328D-02,  1.093D-01/
C
      DATA PDCH /
     *     1.000D+03,  9.453D-02,  9.804D-02,  8.084D-02,  1.208D-01,
     *     6.520D-02,  8.233D-02,  1.084D-01,  1.474D-01,  9.328D-02,
     *     1.093D-01,  1.400D+03,  1.072D-01,  7.450D-02,  6.645D-02,
     *     1.136D-01,  6.750D-02,  8.580D-02,  1.110D-01,  1.530D-01,
     *     1.010D-01,  1.350D-01,  2.170D+03,  4.004D-02,  3.013D-02,
     *     2.664D-02,  5.511D-02,  4.240D-02,  7.660D-02,  1.364D-01,
     *     2.300D-01,  1.670D-01,  2.010D-01,  2.900D+03,  1.870D-02,
     *     1.804D-02,  1.320D-02,  2.970D-02,  2.860D-02,  5.160D-02,
     *     1.020D-01,  2.400D-01,  2.250D-01,  3.370D-01,  4.400D+03,
     *     1.196D-03,  8.784D-03,  1.517D-02,  2.874D-02,  2.488D-02,
     *     4.464D-02,  8.330D-02,  2.008D-01,  2.360D-01,  3.567D-01/
C
      DATA (DCHN(I),I=1,90) /
     *     4.770D-01,  4.750D-01,  4.715D-01,  4.685D-01,  4.650D-01,
     *     4.610D-01,  4.570D-01,  4.550D-01,  4.500D-01,  4.450D-01,
     *     4.405D-01,  4.350D-01,  4.300D-01,  4.250D-01,  4.200D-01,
     *     4.130D-01,  4.060D-01,  4.000D-01,  3.915D-01,  3.840D-01,
     *     3.760D-01,  3.675D-01,  3.580D-01,  3.500D-01,  3.400D-01,
     *     3.300D-01,  3.200D-01,  3.100D-01,  3.000D-01,  2.900D-01,
     *     2.800D-01,  2.700D-01,  2.600D-01,  2.500D-01,  2.400D-01,
     *     2.315D-01,  2.240D-01,  2.150D-01,  2.060D-01,  2.000D-01,
     *     1.915D-01,  1.850D-01,  1.780D-01,  1.720D-01,  1.660D-01,
     *     1.600D-01,  1.550D-01,  1.500D-01,  1.450D-01,  1.400D-01,
     *     1.360D-01,  1.320D-01,  1.280D-01,  1.250D-01,  1.210D-01,
     *     1.180D-01,  1.150D-01,  1.120D-01,  1.100D-01,  1.070D-01,
     *     1.050D-01,  1.030D-01,  1.010D-01,  9.900D-02,  9.700D-02,
     *     9.550D-02,  9.480D-02,  9.400D-02,  9.200D-02,  9.150D-02,
     *     9.100D-02,  9.000D-02,  8.990D-02,  8.900D-02,  8.850D-02,
     *     8.750D-02,  8.700D-02,  8.650D-02,  8.550D-02,  8.500D-02,
     *     8.499D-02,  8.450D-02,  8.350D-02,  8.300D-02,  8.250D-02,
     *     8.150D-02,  8.100D-02,  8.030D-02,  8.000D-02,  7.990D-02/
      DATA (DCHN(I),I=91,143) /
     *     7.980D-02,  7.950D-02,  7.900D-02,  7.860D-02,  7.800D-02,
     *     7.750D-02,  7.650D-02,  7.620D-02,  7.600D-02,  7.550D-02,
     *     7.530D-02,  7.500D-02,  7.499D-02,  7.498D-02,  7.480D-02,
     *     7.450D-02,  7.400D-02,  7.350D-02,  7.300D-02,  7.250D-02,
     *     7.230D-02,  7.200D-02,  7.100D-02,  7.050D-02,  7.020D-02,
     *     7.000D-02,  6.999D-02,  6.995D-02,  6.993D-02,  6.991D-02,
     *     6.990D-02,  6.870D-02,  6.850D-02,  6.800D-02,  6.780D-02,
     *     6.750D-02,  6.700D-02,  6.650D-02,  6.630D-02,  6.600D-02,
     *     6.550D-02,  6.525D-02,  6.510D-02,  6.500D-02,  6.499D-02,
     *     6.498D-02,  6.496D-02,  6.494D-02,  6.493D-02,  6.490D-02,
     *     6.488D-02,  6.485D-02,  6.480D-02/
C
      DATA DCHNA /
     *     6.300D+02,  7.810D-02,  1.421D-01,  1.979D-01,  2.479D-01,
     *     3.360D-01,  5.400D-01,  7.236D-01,  1.000D+00,  1.540D+03,
     *     2.225D-01,  3.950D-01,  5.279D-01,  6.298D-01,  7.718D-01,
     *     9.405D-01,  9.835D-01,  1.000D+00,  2.560D+03,  2.625D-01,
     *     4.550D-01,  5.963D-01,  7.020D-01,  8.380D-01,  9.603D-01,
     *     9.903D-01,  1.000D+00,  3.520D+03,  4.250D-01,  6.875D-01,
     *     8.363D-01,  9.163D-01,  9.828D-01,  1.000D+00,  1.000D+00,
     *     1.000D+00/
C
      DATA DCHNB /
     *     6.300D+02,  3.800D-02,  7.164D-02,  1.275D-01,  2.171D-01,
     *     3.227D-01,  4.091D-01,  5.051D-01,  6.061D-01,  7.074D-01,
     *     8.434D-01,  1.000D+00,  2.040D+03,  1.200D-01,  2.115D-01,
     *     3.395D-01,  5.295D-01,  7.251D-01,  8.511D-01,  9.487D-01,
     *     9.987D-01,  1.000D+00,  1.000D+00,  1.000D+00,  2.200D+03,
     *     1.344D-01,  2.324D-01,  3.754D-01,  5.674D-01,  7.624D-01,
     *     8.896D-01,  9.808D-01,  1.000D+00,  1.000D+00,  1.000D+00,
     *     1.000D+00,  2.850D+03,  2.330D-01,  4.130D-01,  6.610D-01,
     *     9.010D-01,  9.970D-01,  1.000D+00,  1.000D+00,  1.000D+00,
     *     1.000D+00,  1.000D+00,  1.000D+00,  3.500D+03,  3.300D-01,
     *     5.450D-01,  7.950D-01,  1.000D+00,  1.000D+00,  1.000D+00,
     *     1.000D+00,  1.000D+00,  1.000D+00,  1.000D+00,  1.000D+00/
C
C---------------------------------------------------------------
        CST=1.0D+0
C
C*      IS THE KINETIC ENERGY GREATER THAN LIMIT ?
C
        IF (EKIN .GT. 3.5D0) RETURN
C
        IF(KPROJ.EQ.8) GOTO 101
        IF(KPROJ.EQ.1) GOTO 102
C*                                INVALID REACTION
        RETURN
C-------------------------------- NP ELASTIC SCATTERING----------
101   CONTINUE
      IF (EKIN.GT.0.740D+0)GOTO 1000
      IF (EKIN.LT.0.300D+0)THEN
C                                 EKIN .LT. 300 MEV
         IDAT=1
      ELSE
C                                 300 MEV < EKIN < 740 MEV
         IDAT=6
      END IF
C
      ENER=EKIN
      IE=ABS(ENER/0.020D+0)
      UNIV=(ENER-IE*0.020D+0)/0.020D+0
C                                            FORWARD/BACKWARD DECISION
      K=IDAT+5*IE
      BWFW=(DCLIN(K+5)-DCLIN(K))*UNIV + DCLIN(K)
      CALL GRNDM(RNDM,1)
      RND=RNDM(1)
      IF (RND.LT.BWFW)THEN
         VALUE2=-1.0D+0
         K=K+1
      ELSE
         VALUE2=1.0D+0
         K=K+3
      END IF
C
      COEF=(DCLIN(K+5)-DCLIN(K))*UNIV + DCLIN(K)
      CALL GRNDM(RNDM,1)
      RND=RNDM(1)
C
      IF(RND.LT.COEF)THEN
         CALL GRNDM(RNDM,1)
         CST=RNDM(1)
         CST=CST*VALUE2
      ELSE
         CALL GRNDM(RNDM,4)
         R1=RNDM(1)
         R2=RNDM(2)
         R3=RNDM(3)
         R4=RNDM(4)
C
         IF(VALUE2.GT.0.D0)THEN
            CST=MAX(R1,R2,R3,R4)
            GOTO 1500
         ELSE
            CALL GRNDM(RNDM,1)
            R5=RNDM(1)
C
            IF (IDAT.EQ.1)THEN
               CST=-MAX(R1,R2,R3,R4,R5)
            ELSE
               CALL GRNDM(RNDM,2)
               R6=RNDM(1)
               R7=RNDM(2)
               CST=-MAX(R1,R2,R3,R4,R5,R6,R7)
            END IF
C
         END IF
C
      END IF
C
      GOTO 1500
C
C********                                EKIN  .GT.  0.74 GEV
C
1000  CONTINUE
      ENER=EKIN - 0.66D0
C     IE=ABS(ENER/0.02D0)
      IE=ENER/0.02D0
      EMEV=EKIN*1.D+03
C
      UNIV=(ENER-IE*0.020D+0)/0.020D+0
      K=IE
      BWFW=(DCHN(K+1)-DCHN(K))*UNIV + DCHN(K)
      CALL GRNDM(RNDM,1)
      RND=RNDM(1)
C                                        FORWARD NEUTRON
      IF (RND.GE.BWFW)THEN
         DO 1200 K=10,36,9
           IF (DCHNA(K).GT.EMEV) THEN
              UNIVE=(EMEV-DCHNA(K-9))/(DCHNA(K)-DCHNA(K-9))
              CALL GRNDM(RNDM,1)
              UNIV=RNDM(1)
              DO 1100 I=1,8
                 II=K+I
                 P=(DCHNA(II)-DCHNA(II-9))*UNIVE + DCHNA(II-9)
C
                 IF (P.GT.UNIV)THEN
                    CALL GRNDM(RNDM,1)
                    UNIV=RNDM(1)
                    FLTI=I-UNIV
                    GOTO(290,290,290,290,330,340,350,360) I
                 END IF
 1100         CONTINUE
           END IF
 1200    CONTINUE
C
      ELSE
C                                        BACKWARD NEUTRON
         DO 1400 K=13,60,12
            IF (DCHNB(K).GT.EMEV) THEN
               UNIVE=(EMEV-DCHNB(K-12))/(DCHNB(K)-DCHNB(K-12))
               CALL GRNDM(RNDM,1)
               UNIV=RNDM(1)
               DO 1300 I=1,11
                 II=K+I
                 P=(DCHNB(II)-DCHNB(II-12))*UNIVE + DCHNB(II-12)
C
                 IF (P.GT.UNIV)THEN
                   CALL GRNDM(RNDM,1)
                   UNIV=RNDM(1)
                   FLTI=I-UNIV
                   GOTO(120,120,140,150,160,160,180,190,200,210,220) I
                 END IF
 1300          CONTINUE
            END IF
 1400    CONTINUE
      END IF
C
120   CST=1.0D-2*FLTI-1.0D+0
      GOTO 1500
140   CST=2.0D-2*UNIV-0.98D+0
      GOTO 1500
150   CST=4.0D-2*UNIV-0.96D+0
      GOTO 1500
160   CST=6.0D-2*FLTI-1.16D+0
      GOTO 1500
180   CST=8.0D-2*UNIV-0.80D+0
      GOTO 1500
190   CST=1.0D-1*UNIV-0.72D+0
      GOTO 1500
200   CST=1.2D-1*UNIV-0.62D+0
      GOTO 1500
210   CST=2.0D-1*UNIV-0.50D+0
      GOTO 1500
220   CST=3.0D-1*(UNIV-1.0D+0)
      GOTO 1500
C
290   CST=1.0D0 - 2.5D-2*FLTI
      GOTO 1500
330   CST=0.85D0 + 0.5D-1*UNIV
      GOTO 1500
340   CST=0.70D0 + 1.5D-1*UNIV
      GOTO 1500
350   CST=0.50D0 + 2.0D-1*UNIV
      GOTO 1500
360   CST=0.50D0*UNIV
C
1500  RETURN
C
C-----------------------------------  PP ELASTIC SCATTERING -------
C
 102  CONTINUE
      EMEV=EKIN*1.D+03
C
      IF (EKIN.LE.0.500D0) THEN
         CALL GRNDM(RNDM,1)
         RND=RNDM(1)
         CST=2.0D0*RND-1.0D0
         RETURN
C
      ELSE IF (EKIN.LT.1.D0) THEN
         DO 2200 K=13,60,12
            IF (PDCI(K).GT.EMEV) THEN
               UNIVE=(EMEV-PDCI(K-12))/(PDCI(K)-PDCI(K-12))
               CALL GRNDM(RNDM,1)
               UNIV=RNDM(1)
               SUM=0.0D0
               DO 2100 I=1,11
                 II=K+I
                 SUM=SUM + (PDCI(II)-PDCI(II-12))*UNIVE + PDCI(II-12)
C
                 IF (UNIV.LT.SUM)THEN
                   CALL GRNDM(RNDM,1)
                   UNIV=RNDM(1)
                   FLTI=I-UNIV
                   GOTO(55,55,55,60,60,65,65,65,65,70,70) I
                 END IF
 2100          CONTINUE
            END IF
 2200    CONTINUE
      ELSE
         DO 2400 K=12,55,11
            IF (PDCH(K).GT.EMEV) THEN
              UNIVE=(EMEV-PDCH(K-11))/(PDCH(K)-PDCH(K-11))
              CALL GRNDM(RNDM,1)
              UNIV=RNDM(1)
              SUM=0.D0
              DO 2300 I=1,10
                II=K+I
                SUM=SUM + (PDCH(II)-PDCH(II-11))*UNIVE + PDCH(II-11)
C
                IF (UNIV.LT.SUM)THEN
                  CALL GRNDM(RNDM,1)
                  UNIV=RNDM(1)
                  FLTI=UNIV+I
                  GOTO(50,55,60,60,65,65,65,65,70,70) I
                END IF
 2300         CONTINUE
            END IF
 2400    CONTINUE
      END IF
C
50    CST=0.4D0*UNIV
      GOTO 2500
55    CST=0.2D0*FLTI
      GOTO 2500
60    CST=0.3D0 + 0.1D0*FLTI
      GOTO 2500
65    CST=0.6D0 + 0.04D0*FLTI
      GOTO 2500
70    CST=0.78D0 + 0.02D0*FLTI
C
2500  CONTINUE
      CALL GRNDM(RNDM,1)
      RND=RNDM(1)
      IF (RND.GT.0.5D+00)CST=-CST
C
      RETURN
      END
